	%%%%%%%%%%%%%%%%%%%%%%%%%%%
	% run3_master.m
	%%%%%%%%%%%%%%%%%%%%%%%%%%%

clear
close all

    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % SET DIRECTORY
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    
	% specify directory:
cd code

output_directory='./../output/';
input_directory='./../input/';
results_directory='./../results/';

MainRun=1;
ExogTest=1; 
DWLboot=1; 

	%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
	% INITIALIZE 
	%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 
fprintf(['\n -------------' datestr(clock) '------------- \n']);

p = gcp('nocreate'); 
if isempty(p)
	parpool(32);
end

rng(123456,'twister')

	%%%%%%%%%%%%%%%%%%%%%%%%%%%
	% CALL ESTIMATION 
	%%%%%%%%%%%%%%%%%%%%%%%%%%%

if MainRun

		% MAIN RUN

	berk_ml.main_compute('state',0,3,7,3,1,0,input_directory,output_directory); 
	
		% INCREASE PRICE POLYNOMIAL

	berk_ml.main_compute('state',0,4,7,3,1,1,input_directory,output_directory); 

		% INCREASE Q POLYNOMIAL
	
	berk_ml.main_compute('state',0,3,8,3,1,0,input_directory,output_directory); 

		% INCREASE Y POLYNOMIAL

	berk_ml.main_compute('state',0,3,7,4,1,0,input_directory,output_directory); 

		% VARY BERKSON ERROR

	berk_ml.main_compute('state',0,3,7,3,2,0,input_directory,output_directory); 
	berk_ml.main_compute('state',0,3,7,3,3,0,input_directory,output_directory); 
	
		% USE EMPIRICAL CDF

	berk_ml.main_compute('state',1,3,7,3,1,0,input_directory,output_directory); 

	fprintf('Note: computational step complete.\n');

		% MAIN RUN

	berk_ml.main_output('state',0,3,7,3,1,output_directory,1,1,1,1,0); 

		% VARY BERKSON
	
	berk_ml.main_output('state',0,3,7,3,2,output_directory,1,0,0,0,0); 
	berk_ml.main_output('state',0,3,7,3,3,output_directory,1,0,0,0,0); 

		% INCREASE PRICE POLYNOMIAL

	berk_ml.main_output('state',0,4,7,3,1,output_directory,0,0,0,1,1); 

		% INCREASE Q POLYNOMIAL
	
	berk_ml.main_output('state',0,3,8,3,1,output_directory,0,0,0,1,1); 
	
		% INCREASE Y POLYNOMIAL

	berk_ml.main_output('state',0,3,7,4,1,output_directory,0,0,0,1,1); 

		% USING STATE-VARYING EMP CDF

	berk_ml.main_output('state',1,3,7,3,1,output_directory,0,0,0,1,1); 

end % MainRun

if ExogTest

		%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
		% START EXOGENEITY TEST 
		%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
	
	K1=3;
	K2=7;
	trim_perc_q=1;

	sigma_berkson_base=0.033;	
	sigma_berkson_high=1.2*0.033;	
	sigma_berkson_low=0.8*0.033;	

	fname=[output_directory 'table-4.txt'];
	fileID=fopen(fname,'w');

	run_exog_test_NHTS(65000,80000,sigma_berkson_base,K1,K2,trim_perc_q,input_directory,output_directory,fileID);
	run_exog_test_NHTS(65000,80000,sigma_berkson_low,K1,K2,trim_perc_q,input_directory,output_directory,fileID);
	run_exog_test_NHTS(65000,80000,sigma_berkson_high,K1,K2,trim_perc_q,input_directory,output_directory,fileID);

	run_exog_test_NHTS(50000,65000,sigma_berkson_base,K1,K2,trim_perc_q,input_directory,output_directory,fileID);
	run_exog_test_NHTS(50000,65000,sigma_berkson_low,K1,K2,trim_perc_q,input_directory,output_directory,fileID);
	run_exog_test_NHTS(50000,65000,sigma_berkson_high,K1,K2,trim_perc_q,input_directory,output_directory,fileID);
	
	run_exog_test_NHTS(35000,50000,sigma_berkson_base,K1,K2,trim_perc_q,input_directory,output_directory,fileID);
	run_exog_test_NHTS(35000,50000,sigma_berkson_low,K1,K2,trim_perc_q,input_directory,output_directory,fileID);
	run_exog_test_NHTS(35000,50000,sigma_berkson_high,K1,K2,trim_perc_q,input_directory,output_directory,fileID);

	fclose(fileID);

end % ExogTest

if DWLboot
	
	constr=0;
	mat_type='chebyshev'; 
	trim_perc_q=1;
	y_min=1;
	y_max=999000;
	range_extend=0.033;
	K1=3;
	K2=7;
	K3=3;
	RR=499;   
	factor_berkson=1;
	Berkson_estimated=0;
	berkson_ref='state';
	berkson_txt=[berkson_ref '_hist' num2str(Berkson_estimated)];

		%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
		% COMPUTE DWL
		%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
	 
	fname=[output_directory 'table-A1.txt'];
	fileID=fopen(fname,'w');

	fprintf(fileID,'\n ---------------------------------\n WITH BERKSON ERRORS, BERKSON TYPE %20s:\n ---------------------------------\n',berkson_txt);
	compute_DWL_bootstrap(mat_type,range_extend,K1,K2,K3,trim_perc_q,y_min,y_max,constr,RR,berkson_txt,factor_berkson,0,input_directory,output_directory,fileID); 

	fprintf(fileID,'\n ---------------------------------\n WITHOUT BERKSON ERRORS, BERKSON TYPE %5s:\n ---------------------------------\n',berkson_txt);
	compute_DWL_bootstrap(mat_type,range_extend,K1,K2,K3,trim_perc_q,y_min,y_max,constr,RR,berkson_txt,factor_berkson,1,input_directory,output_directory,fileID); 


end % DWLboot

	%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
	% COPY FILES
	%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 
copyfile([output_directory 'fig_tau50_ff1_yy57500_BW.pdf'],[results_directory 'figure-4.pdf']);
copyfile([output_directory 'fig_combined_ff1_yy57500_BW.pdf'],[results_directory 'figure-5.pdf']);
copyfile([output_directory 'legend_BW.pdf'],[results_directory 'figure-6_legend.pdf']);
copyfile([output_directory 'fig_tau50_ff1_yy72500_slim_BW.pdf'],[results_directory 'figure-6-a.pdf']);
copyfile([output_directory 'fig_tau50_ff1_yy57500_slim_BW.pdf'],[results_directory 'figure-6-b.pdf']);
copyfile([output_directory 'fig_tau50_ff1_yy42500_slim_BW.pdf'],[results_directory 'figure-6-c.pdf']);
copyfile([output_directory 'fig_tau50_ff1_yy57500_BW.pdf'],[results_directory 'figure-7-a.pdf']);
copyfile([output_directory 'fig_tau50_ff3_yy57500_BW.pdf'],[results_directory 'figure-7-b.pdf']);
copyfile([output_directory 'fig_tau50_ff2_yy57500_BW.pdf'],[results_directory 'figure-7-c.pdf']);

copyfile([output_directory 'DWL_state_hist0_ff1_allY_p3_q7_y3.txt'],[results_directory 'table-3.txt']);
copyfile([output_directory 'table-4.txt'],[results_directory 'table-4.txt']);
copyfile([output_directory 'table-A1.txt'],[results_directory 'table-A1.txt']);
copyfile([output_directory 'DWL_state_hist0_ff1_allY_p4_q7_y3.txt'],[results_directory 'table-A2-a.txt']);
copyfile([output_directory 'DWL_state_hist0_ff1_allY_p3_q7_y4.txt'],[results_directory 'table-A2-b.txt']);
copyfile([output_directory 'DWL_state_hist0_ff1_allY_p3_q8_y3.txt'],[results_directory 'table-A2-c.txt']);
copyfile([output_directory 'DWL_state_hist0_ff1_allY_p3_q7_y3.txt'],[results_directory 'table-A3-a.txt']);
copyfile([output_directory 'DWL_state_hist1_ff1_allY_p3_q7_y3.txt'],[results_directory 'table-A3-b.txt']);

	%%%%%%%%%%%%%%%%%%%%%%%%%%%
	% END OF PROGRAM 
	%%%%%%%%%%%%%%%%%%%%%%%%%%%

p = gcp('nocreate'); 
if isempty(p)==0
	delete(p)
end

fprintf(['\n -------------' datestr(clock) '------------- \n']);

quit

	%
